library(sf)
library(tidyverse)
library(stringi)
library(DT)
library(ggplot2)
library(tmap)Script projet CIV6708 H23
Script Quarto ayant pour but de documenter les traitements effectués dans le cadre du projet CIV6708.
Date du dernier render : 2023-03-24
1 Librairies
2 Recensement
2.1 Données
2.1.1 Population et tranches d’âge
Données brut du recensement avec la population par tranches d’âge, par Aires de diffusion (AD).
Source : recensement 2021
recensement.age <- read.csv("C:/Users/victo/OneDrive - polymtl.ca/CIV6708/PROJET/StatCan/recensement_age_AD.csv",
sep = ";", fileEncoding = "UTF-8-BOM")Echelle : Canada entier (907752 lignes).
Affichage des premières lignes :
2.1.2 Limite des AD
On charge les limites géographiques des AD de Drummondville :
sf.drummond.ville.AD <- st_read("C:/Users/victo/OneDrive - polymtl.ca/CIV6708/PROJET/SDR_drummondville/AD_drummondville.shp")Reading layer `AD_drummondville' from data source
`C:\Users\victo\OneDrive - polymtl.ca\CIV6708\PROJET\SDR_drummondville\AD_drummondville.shp'
using driver `ESRI Shapefile'
Simple feature collection with 137 features and 4 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 7684020 ymin: 1299056 xmax: 7719001 ymax: 1324995
Projected CRS: NAD83 / Statistics Canada LambertAffichage des éléments :
2.2 Traitement
On commence par formater les noms des colonnes :
colnames <- colnames(recensement.age)
colnames <- stri_trans_general(colnames, id = "Latin-ASCII") # retire les accents
colnames <- tolower(colnames) # met tout en minuscules
colnames <- str_replace_all(colnames, "\\.", "_") # remplace les points par des underscores
colnames(recensement.age) <- colnamesOn transforme les données :
recensement.age <- recensement.age %>%
select( # sélection des colonnes intéressantes
c(
geo,
dguid,
age__16_,
valeur
)
) %>%
filter( # sélection des données de Drummondville seulement
geo %in% sf.drummond.ville.AD$ADIDU
) %>%
pivot_wider( # on pivote les données dans la largeur (afin d'avoir une ligne par AD)
names_from = "age__16_", values_from = "valeur"
)On procède aux mêmes traitements que précédemment sur les nouvelles colonnes :
colnames <- colnames(recensement.age)
colnames <- stri_trans_general(colnames, id = "Latin-ASCII")
colnames <- tolower(colnames)
colnames <- str_replace_all(colnames, "\\.", "_")
colnames <- str_replace_all(colnames, " ", "_")
colnames <- str_remove(colnames, "-")
colnames(recensement.age) <- colnamesOn remplace les NA par des 0 :
recensement.age[is.na(recensement.age)] <- 0On affiche les données traitées :
2.3 Enregistrement
write.csv(recensement.age, "C:/Users/victo/OneDrive - polymtl.ca/CIV6708/PROJET/StatCan/recensement_age_AD_tri.csv", row.names = F)3 Association de la population aux bâtiments
3.1 Données
On charge les données de population par tranches d’âge (traitées précédemment) :
recensement.age <- read.csv("C:/Users/victo/OneDrive - polymtl.ca/CIV6708/PROJET/StatCan/recensement_age_AD_tri.csv")On charge également les bâtiments résidentiels, extraits de OSM.
Clé utilisée : building
Valeurs : apartments, cabin, detached, dormitory, semidetached_house, static_caravan
sf.accomodations <- st_read("C:/Users/victo/OneDrive - polymtl.ca/CIV6708/PROJET/OSM/buildings_accomodations_drummondville.shp")Reading layer `buildings_accomodations_drummondville' from data source
`C:\Users\victo\OneDrive - polymtl.ca\CIV6708\PROJET\OSM\buildings_accomodations_drummondville.shp'
using driver `ESRI Shapefile'
Simple feature collection with 5329 features and 31 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -72.64695 ymin: 45.75518 xmax: -72.28676 ymax: 46.02368
Geodetic CRS: WGS 84Avant d’afficher les données, on fait quelques sélections et changements dans les noms des champs :
sf.accomodations <- sf.accomodations %>%
select( # sélection des champs qui nous intéressent
c(
full_id:building,
building_1:geometry
)
) %>%
rename(
flats = building_1 # on corrige l'abbréviation automatique du format shapefile
)Affichage :
Pour information, le champ acc_number a été créé dans QGIS, et vaut 1 si flats est vide, 0 sinon.
3.2 Traitement
On commence par changer le CRS des données OSM, pour concorder avec celui des données de StatCan :
sf.accomodations <- st_transform(sf.accomodations, crs = 3347)On associe à chaque bâtiment le code de l’AD dans laquelle son centroïde se trouve :
sf.accomodations.points <- st_join(st_centroid(sf.accomodations), sf.drummond.ville.AD[c("IDUGD")],
join = st_within)Warning: st_centroid assumes attributes are constant over geometriessf.accomodations$IDUGD <-
sf.accomodations.points$IDUGD[match(sf.accomodations$full_id,
sf.accomodations.points$full_id)]On ajoute à chaque bâtiment les données de population pour l’AD dans laquelle il se trouve :
sf.accomodations <- left_join(sf.accomodations, recensement.age,
by = c("IDUGD" = "dguid"))Affichage des premières lignes de la table des bâtiments :
On a actuellement, pour chaque bâtiment, les données de population pour l’AD au complet. On souhaite répartir ces données à l’échelle de chaque bâtiment, en se basant sur deux conditions :
- les sommes des données de population pour tous les bâtiments d’une AD doit être égale aux données de StatCan ;
- le nombre de personnes attribuées à chaque bâtiment doit être proportionnel au nombre de logements de ce bâtiments.
On va pour cela utiliser la fonction suivante, qui permet de générer un vecteur de nombre entier à partir de la somme de ce vecteur et d’un vecteur de probabilités.
rand_vect <- function(N, M, sd = 1, pos.only = TRUE,prob) {
vec <- rbinom(N,1,prob)
deviation <- M - sum(vec)
for (. in seq_len(abs(deviation))) {
vec[i] <- vec[i <- sample(N, 1,replace=TRUE,prob)] + sign(deviation)
}
if (pos.only) while (any(vec < 0)) {
negs <- vec < 0
pos <- vec > 0
vec[negs][i] <- vec[negs][i <- sample(sum(negs), 1)] + 1
vec[pos][i] <- vec[pos ][i <- sample(sum(pos ), 1)] - 1
}
vec
}On attribue les personnes à chaque bâtiment :
df.accomodations <- st_drop_geometry(sf.accomodations) %>% # on travaille avec un dataframe plutôt qu'un sf
group_by(IDUGD) %>% # on travaille par AD
mutate(
nb_bat_AD = n(), # nombre de bâtiments dans l'AD
share_log_AD = round(acc_number/sum(acc_number), 3), # part des logements de l'AD représentée par le bâtiment
total_pop_AD = total__age,
across(X15_a_19_ans:X85_ans_et_plus, # opération effectuée à l'identique pour chaque tranche d'âge
~ rand_vect(N = n(), # nombre de bâtiments à populer dans l'AD
M = unique(.x), # nombre de personne dans l'AD pour cette tranche d'âge
prob = share_log_AD) # proportion des logements représentée par chaque bâtiment pour l'AD
)
) %>%
rowwise() %>% mutate(
total__age = sum(c_across(X15_a_19_ans:X85_ans_et_plus)) # pour le total, on somme ce qu'on a obtenu pour chaque tranche d'âges
) %>%
# quelques opérations pour réorganiser les colonnes :
select(-geo) %>%
relocate(IDUGD) %>%
relocate(nb_bat_AD,
share_log_AD,
total_pop_AD,
.before = total__age)
sf.accomodations <- left_join(sf.accomodations[c("full_id")], df.accomodations) %>%
relocate(IDUGD) # on ajoute les traitements aux géométries d'origineJoining with `by = join_by(full_id)`rm(df.accomodations)Affichage de la table complète :
Warning in instance$preRenderHook(instance): It seems your data is too big for
client-side DataTables. You may consider server-side processing:
https://rstudio.github.io/DT/server.htmlCarte :
La suite des traitements (production d’une carte de chaleur) va être effectuée dans QGIS.
3.3 Enregistrement
st_write(sf.accomodations, "C:/Users/victo/OneDrive - polymtl.ca/CIV6708/PROJET/OSM/accomodations_drummondville.shp",
delete_layer = T)Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
Shapefile driverDeleting layer `accomodations_drummondville' using driver `ESRI Shapefile'
Writing layer `accomodations_drummondville' to data source
`C:/Users/victo/OneDrive - polymtl.ca/CIV6708/PROJET/OSM/accomodations_drummondville.shp' using driver `ESRI Shapefile'
Writing 5329 features with 26 fields and geometry type Multi Polygon.